# Code to create Figure A5 in "Pandemics and Cities: Evidence from the Black 
# Death and the Long-Run"
# Author: Noel Johnson
# Date Created: 12-10-21
# Last updated: 9-8-23

library(tidyverse)
library(sf)
library(tmap)
library(haven)

# set working directory
setwd("/Users/noeljohnson_laptop/Dropbox/Research/Plague City Growth/JUE RR1 Nov 2022/figures_replication/Figure_A5")

# read in, project, and plot the 1300 borders map from Nussli
map_1300 <- st_read("gis_1300/data/1300/latin1/sovereign_dependent_states.shp")
map_1300 <- st_transform(map_1300, crs=4326)
map_1300
glimpse(map_1300)
plot(map_1300[12], axes=T, reset=T)

# read in, project, and plot the cities
cities_165 <- read_csv("cities_for_map165.csv")
# spatialize the city data. turn into sf points, specify coords and crs
cities_165_sf <- st_as_sf(cities_165,
                      coords = c("longitude", "latitude"),
                      crs = 4326)
plot(cities_165_sf[4], axes=T, reset=T)

# plot the borders and cities together
plot(map_1300[12], reset=F)
plot(cities_165_sf[4], add=T)

# need to crop the map and make it look prettier 

# Make a bounding box of the cities so you can create one a little bigger
bbox_cities <- st_bbox(cities_165_sf, crs=4326)
bbox_cities

xrange <- bbox_cities$xmax - bbox_cities$xmin # range of x values
yrange <- bbox_cities$ymax - bbox_cities$ymin # range of y values

bbox_cities[1] <- bbox_cities[1] - (0.1 * xrange) # xmin - left
bbox_cities[3] <- bbox_cities[3] + (0.1 * xrange) # xmax - right
bbox_cities[2] <- bbox_cities[2] - (0.1 * yrange) # ymin - bottom
bbox_cities[4] <- bbox_cities[4] + (0.1 * yrange) # ymax - top

bbox_cities_sf <- bbox_cities %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

bbox_cities_sf
# plot(bbox_cities_sf)

# plot the borders, cities, and bounding box together
plot(map_1300[12], reset=F)
plot(cities_165_sf[4], add=T)
plot(bbox_cities_sf, add=T)

# clip the map and cities using the bounding box
cities_165_crop <- st_crop(cities_165_sf, bbox_cities_sf)

# have to fix the geometry of map_1300 before it can be cropped
library(lwgeom)
st_is_valid(map_1300)
map_1300_fixed <- st_make_valid(map_1300)
st_is_valid(map_1300_fixed) # this doesn't work
sf::sf_use_s2(FALSE) # this does work
map_1300_crop <- st_crop(map_1300, bbox_cities_sf)
sf::sf_use_s2(TRUE) # putting settings back to their default

# plot the cropped borders and cities
plot(map_1300_crop[12], reset=F)
plot(cities_165_crop[4], add=T)

# Make the figure for the paper...
state_map <- tm_shape(map_1300_crop) + tm_borders(col = "black") 
state_map

mycols <- colors()[c(142, 54, 35, 41, 24)]
cities_map_color <- tm_shape(cities_165_crop) +
  tm_bubbles("mortality", col = "mortality", border.col= "black",
             border.alpha = 0.5, style="fixed",
             breaks = c(-Inf, 20, 40, 60, 80, Inf),
             palette=mycols, title.size = "Mortality Rate (Size)",
             scale = 2.0, title.col ="Mortality Rate (Color)") +
  tm_layout(legend.stack = "vertical", legend.text.size = .80,
            frame = F, inner.margins = .10,
            legend.position = c(-.15, 0.05)) +
  state_map

cities_map_color

# save map
tmap_save(cities_map_color, filename="cities_map_color.png",
          height=8.5, width=11, units="in", dpi=300)

#







# End Code







